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ABSTRACT 

We model the variability profiles of millisecond period X-ray pulsars. We performed three- 
dimensional magnetohydrodynamic simulations of disk accretion to millisecond period neutron stars 
with a misaligned magnetic dipole moment, using the pseudo-Newtonian Paczyhski-Wiita potential to 
model general relativistic effects. We found that the shapes of the resulting funnel streams of accreting 
matter and the hot spots on the surface of the star are quite similar to those for more slowly rotating 
stars obtained from earlier simulations using the Newtonian potential. The funnel streams and hot 
spots rotate approximately with the same angular velocity as the star. The spots are bow-shaped 
(bar-shaped) for small (large) misalignment angles. We found that the matter falling on the star has 
a higher Mach number when we use the Paczyhski-Wiita potential than in the Newtonian case. 

Having obtained the surface distribution of the emitted flux, we calculated the variability curves 
of the star, taking into account general relativistic, Doppler and light-travel-time effects. We found 
that general relativistic effects decrease the pulse fraction (flatten the light curve), while Doppler and 
light-travel-time effects increase it and distort the light curve. We also found that the light curves 
from our hot spots are reproduced reasonably well by spots with a gaussian flux distribution centered 
at the magnetic poles. We also calculated the observed image of the star in a few cases, and saw that 
for certain orientations, both the antipodal hot spots are simultaneously visible, as noted by earlier 
authors. 

Subject headings: accretion, accretion disks — pulsars — X-rays: stars 



1. INTRODUCTION 

Millisecond X-ray pulsars show bursts of periodic and 
quasi-periodic variability in the X-ray (Stella & Vietri 
1999; van der Klis 2000; Chakrabarty et al. 2003; Wij- 
nands et al. 2003). Six of these pulsars are believed to 
be accretion-powered (Wijnands 2005). Numerical mod- 
elling of accretion is an important tool for studying these 
phenomena. Three-dimensional magnetohydrodynamic 
(3D MHD) simulations of disk accretion to rotating mag- 
netized stars with a misaligned dipole magnetic field have 
been carried out by Koldoba et al. (2002) and Romanova 
et al. (2003; 2004, hereafter - R04). They found that 
the accreting matter is channelled by the star's magnetic 
field into two antipodal streams or funnels, and falls on 
the stellar surface forming two antipodal hot spots - re- 
gions of relatively high temperature. Such flows have 
been predicted theoretically (e.g., Pringle & Rees 1972; 
Lamb, Pethick & Pines 1973; Ghosh & Lamb 1978, 1979), 
but were modelled numerically only recently. Assuming 
that the energy of the infalling matter is converted en- 
tirely into radiation, R04 then calculated the variability 
curves of the star. Those simulations and calculations 
were performed for a generic star in a completely clas- 
sical framework. We refined those simulations, focusing 
on rapidly rotating (~ 3-5 ms period) neutron stars. For 
such stars, the following issues arise: (1) For compact 
objects like neutron stars, general relativistic effects sig- 
nificantly influence the accretion process as well as the 
observed flux, (i) We modelled general relativistic effects 



on accretion by using the pseudo-Newtonian Paczyhski- 
Wiita potential (Paczyhski and Wiita 1980) which re- 
produces some important features of the Schwarzschild 
geometry, like the positions of the innermost stable and 
marginally bound circular orbits, (ii) The variability 
curve of the star also changes significantly, because grav- 
itational bending of light emitted by the star allows more 
of the star than the hemisphere facing the observer to be 
visible, and gravitational redshift of the light decreases 
the total flux observed. We use the Schwarzschild metric 
to take these effects into account. (2) For rapidly rotat- 
ing stars like millisecond pulsars, the rotation of the star 
changes the observed flux through the twin special rela- 
tivistic effects of Doppler shift and relativistic beaming of 
the emitted radiation. Henceforth we refer to these two 
effects collectively as the "Doppler effect." (3) The time 
difference between light emitted from different points on 
the star reaching the observer is important when the lin- 
ear speed of the emitting region is comparable to the 
speed of light, and also when the emitting object is com- 
pact, and causes distortion of the observed shape of the 
hot spots (the apparent position of a point on the stel- 
lar surface can differ by as much as 10° from its actual 
position) . We build upon the earlier calculations to take 
these effects into account. Additionally, due to the high 
rotation speed, the Kerr metric would be closer to the 
actual metric around the star than the Schwarzschild 
metric. However we do not use the Kerr metric, since 
we find from numerical integrations that the frame drag- 
ging effects introduced by the Kerr metric are relatively 
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small (see also Braje, Romani & Rauch 2000), particu- 
larly when compared with the other errors introduced by 
the assumptions in our variability model. 

These effects on the light curves have been taken into 
account by earlier authors to obtain light curves for sim- 
ple hot spots (see, e.g., Pechenick, Ftaclas & Cohen 1983; 
Ftaclas, Kearney & Pechenick 1986; Braje, Romani & 
Rauch 2000; Ford 2000; Beloborodov 2002; Poutanen 
& Gierlinski 2003; Viironen & Poutanen 2004; Bhat- 
tacharyya et al. 2005). Here, we obtain light curves 
using realistic hot spots obtained from our accretion sim- 
ulations. 

In section 2 we present the results of our MHD simu- 
lations. We then discuss the analytical background for 
calculating variability curves in section 3, and show the 
variability curves for synthetic and realistic spots in sec- 
tions 4 and 5 respectively, followed by some concluding 
remarks in section 6. 

2. DISK ACCRETION - 3D SIMULATIONS 

2.1. Earlier Simulations 

We briefly describe earlier 3D MHD simulations 
(Koldoba et al. 2002; Romanova et al. 2003; R04). The 
star has a dipolc magnetic field, the axis of which makes 
an angle O with the star's rotation axis. The rotation 
axes of the star and the accretion disk are aligned. The 
disk has a low-density corona which also rotates about 
the same axis. To model stationary accretion, the disk 
was chosen to initially be in a quasi-cquilibrium state, 
where the gravitational, centrifugal and pressure gradi- 
ent forces are in balance (Romanova et al., 2002). Vis- 
cosity is modelled using the a-model (Novikov & Thorne 
1973; Shakura & Sunyaev 1973). To model accretion, the 
ideal MHD equations were solved numerically in three 
dimensions, using a Godunov-type numerical code, writ- 
ten in a "cubed-sphere" coordinate system rotating with 
the star (Koldoba et al., 2002; Romanova et al., 2003). 
The boundary conditions at the star's surface amount 
to assuming that the infalling matter passes through the 
surface of the star. So the dynamics of the matter after 
it falls on the star was ignored. It was found that the in- 
ward motion of the accretion disk is stopped by the star's 
magnetosphere at the Alfven radius, where the magnetic 
and matter energy densities become equal. At that point 
the matter leaves the disk and moves along the magnetic 
field lines. This flow is called a funnel stream. In this 
region, the matter radiates primarily in the X-ray. It 
heats up the star's surface where it falls, forming "hot 
spots." There are two antipodal funnel streams and hot 
spots. From the point of view of an external observer, 
they rotate with approximately the same angular veloc- 
ity as the star, causing the observed flux to vary period- 
ically with time. The shape of the funnel streams and 
hot spots keeps changing slightly with time, leading to 
quasi-variability in the observed flux. 

2.2. Reference Values 

In our new simulations, we use the same model as 
described above. The simulations are done using the 
following dimensionless variables: the radial coordinate 
r 1 = r/Ro, the fluid velocity v' = v/vq, the density 
p' = p/po, the magnetic field B' = B/_B , the pressure 
p' = P/Po, the temperature T" = T/Tq, and the time 



t' = t/to. The variables with subscript arc dimensional 
reference values and the unprimed variables are the di- 
mensional variables. Because of the use of dimensionless 
variables, the results are applicable to a wide range of 
objects and physical conditions, each with its own set 
of reference values. To apply our simulation results to 
a particular situation, we have the freedom to choose 
three parameters, and all the reference values are calcu- 
lated from those. We choose the mass, radius and surface 
magnetic field of the star as the three independent pa- 
rameters. 

The reference values are determined as follows: The 
unit of distance Rq is chosen such that the star has radius 
R = 0.35i?o- The reference velocity is the Keplerian 
velocity at Rq, v = {GM / Rq) 1 / 2 , and uj = v /Ro is the 
reference angular velocity. The reference time is to = 
Ro/vq. The reference surface magnetic field of the star 
is B+ . The reference magnetic field, B , is the initial 
magnetic field strength at r = Rq, assuming a surface 
magnetic field of B± a . The reference density is taken to 
be po = B 2 /v 2 . The reference pressure is po = Pov 2 - The 
reference temperature is T = po /Tlpo, where 1Z is the gas 
constant. The reference accretion rate is M = PqVqRq. 
The reference energy flux is Eq = p^v^R 2 ,. The reference 
value for the effective blackbody temperature of the hot 
spots is (T e ff) = (poVq/o-) 1 ^ 4 , where a is the Stefan- 
Boltzmann constant. 

For the millisecond pulsars in our simulations, we take 
the mass of the neutron star to be M = 1.4M Q = 
2.8 x 10 33 g and its radius R = 10 km = 10 6 cm. The 
reference length scale is Rq w 2.86i? = 2.86 x 10 6 cm. 
The reference velocity is vo = 8.1 x 10 9 cm s _1 . The 
reference time is to — 0-35 ms. The reference surface 
magnetic field is _B i0 = 10 8 G, which is a typical value 
for millisecond pulsars. Then the reference magnetic 
field is B = B+ (R/R ) 3 ~ 4.3 x 10 6 G. The ref- 
erence density is p — 2.8 x 10~ 7 g cm~ 3 . The ref- 
erence pressure is po = 1.8 x 10 13 dynes cm~ 2 . The 
reference temperature is To = 7.9 x 10 11 K. The ref- 
erence value of the effective blackbody temperature is 
(T e ff) « 7.2 x 10 6 K. The reference mass accretion rate 
is M « 1-85 x 10 16 g s- 1 w 2.9 x 1O- 1O M yr" 1 . The 
reference energy flux is Eq w 1.2 x 10 36 erg s _1 . 

Subsequently, we drop the primes on the dimensionless 
variables and show dimensionless values in the figures. 

2.3. New Simulations 

We performed new simulations for rapidly rotating (3-5 
ms period) neutron stars. We followed the same proce- 
dure as described above, with two major changes: (1) 
We modelled general relativistic effects on accretion by 
using the pseudo-Newtonian Paczynski-Wiita (PW) po- 
tential (Paczyhski & Wiita 1980), $(r) = -GM/(r-r g ), 
where M is the mass of the star and r g = 2GM/c 2 is its 
Schwarzschild radius. (2) We calculated the variability 
curves taking into account relativistic and light-travel- 
timc effects. 

We used the following parameters in our simulations: 
The surface magnetic field of the star was 5 x 10 7 G. 
The star's Schwarzschild radius was r g = 4.15 km = 
4.15 x 10 5 cm, or 0.145 in dimensionless units. The disk 
and corona initially had temperatures of 0.01 and 1 re- 
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Fig. 1. — Mass accretion rates for a 3-ms pulsar with (solid lines) and without (dotted lines) the Paczynski-Wiita potential, for different 
misalignment angles. 
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Fig. 2. — Top row: Matter flow around a star with P = 3ms in the "PW" case for different misalignment angles. The red lines are 
sample magnetic field lines and the white lines are sample streamlines of matter flow. Second row: Distribution of emitted flux on the star's 
surface, without using the PW potential. Third row: Distribution of emitted flux on the star's surface, using the PW potential. Bottom 
row: Distribution of emitted flux on a sphere of radius 15 km, using the PW potential. 



spectively, and densities of 1 and 0.01 respectively. The 
disk was thinner than in the earlier simulations. The vis- 
cosity a-parameter was 0.04. Each of the six blocks of 
our cubed-sphere grid had 65 cells in the radial direction 



and 31 in the angular direction, which corresponds to an 
outer disk radius of ~ 8.7, or about 25 stellar radii. 

We used a modified code and modified initial condi- 
tions that use the PW potential instead of the Newtonian 
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one. The initial values of physical variables in the disk 
do not change significantly due to use of the PW poten- 
tial. We found that the shapes of the funnel streams and 
hot spots are similar to those in earlier simulations for 
more slowly rotating stars, and in simulations with the 
Newtonian potential. We also found that the velocity of 
the infalling matter near the surface of the star is higher 
in the "PW" case than in the "non-PW" case. This is 
expected because the PW potential is stronger than the 
Newtonian one. 

The accretion rate is higher in the PW case than in the 
non-PW Fig. n shows. This is again expected 

because the PW potential is stronger. This will lead 
to faster depletion of the inner disk matter. If angular 
momentum is efficiently transported outward, e.g., by 
the magneto-rotational instability (see, e.g., Hawley & 
Krolik 2001), then the influence of the PW potential will 
eventually be felt in the outer disk regions, and enhanced 
accretion could be sustained for a longer time. 

The top row of Fig. |2Jshows the flow of matter around 
a star with P = 3 ms in the PW case. The flow is similar 
in the non-PW case. The next two rows show the hot 
spots formed on the surface of the star, without and with 
the PW potential. We see that the position and shape 
of the hot spots depends on the misalignment angle, but 
does not significantly depend on the presence of the PW 
potential. The hot spots are bow shaped (bar-shaped) 
for small (large) misalignment angles. The emitted flux 
is highest at the center of the spots, and decreases out- 
wards. Note that the hot spots are not usually centered 
at the magnetic poles. In fact, they do not even fall 
on the magnetic poles in most cases. We see that the 
emitted flux is higher in the PW case, which is expected 
because the PW potential is stronger than the Newtonian 
one. 

To get an idea of the conditions at the surface of a 
larger neutron star with the same mass, rotation period 
and magnetic dipole moment, in a similar situation, we 
can look at the surface of a sphere of radius 15 km con- 
centric with the star in our simulations. This approach is 
valid because changing the star's radius does not change 
the accretion flow around the star, since the PW po- 
tential depends only on the star's Schwarzschild radius. 
The bottom row of Fig. El shows the hot spots on such 
a sphere. We see that, other conditions remaining the 
same, a larger star has much fainter hot spots, which is 
again to be expected because the accreting matter has a 
lower velocity at the surface of the larger star. 

For 3-ms pulsars, the typical values of physical quan- 
tities observed in our simulations after 4-6 rotations of 
the star are as follows: The surface magnetic field does 
not change in our model, and hence is ~ 5 x 10 7 G. The 

mass accretion rate to the star ~ 10 15 g s^ 1 ~ 10 _12 M, 
-l 



(a) 





yr~". The total power emitted from the star, after cor- 
recting for gravitational redshift, ~ 10 34 erg s _1 . In the 
hot spots, the matter density ~ 10~ 7 g cm -3 . The speed 
of the inflowing matter ~ 2 xlO 10 cm s _1 . The effec- 
tive blackbody temperature of the hot spots ~ 5 x 10 6 
K. The matter pressure ~ 10 12 dynes cm -2 . The Mach 
number ~ 1 — 6 in the PW case and ~ 0.8 — 3 in the 
non-PW case. The surface distributions of the density, 
pressure, velocity and temperature closely follow that of 
the emitted flux. 
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Fig. 3. — Path of light emitted by a compact object: (a) neglect- 
ing light bending, and (b) including light bending. 



3. CALCULATION OF THE VARIABILITY CURVES 

Having obtained the distribution of emitted flux on the 
star's surface, we can calculate the flux received by an 
observer. Without including the relativistic and light- 
travel-time effects, the observed flux at a large distance 
D from the star is proportional to 



J = 



dS Ib(R, ip) cosip, 



(1) 



cos i/>>0 



where dS is an element of area of the star's surface at 
a position R, and Ie (R, VO is t ne intensity of radiation 
emitted by dS at an angle tp with respect to the local 
radial direction (Fig. |2t). We have ignored the overall 
factor of 1/ D 2 in the flux. When calculating this numeri- 
cally, the integral becomes a sum over grid elements. We 
now consider the gravitational and rotational effects one 
by one. These calculations have been done by earlier au- 
thors (see, e.g., Poutanen & Gierlinski 2003; Viironen & 
Poutanen 2004), but we present them here for the sake 
of completeness and for discussing slight differences in 
approach. 

3.1. General Relativistic Effects 

When taking general relativistic effects into account, 
two approaches are possible towards calculating the ob- 
served flux. In the traditional ray tracing method (see, 
e.g., Braje, Romani & Rauch 2000; Bhattacharyya et 
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al. 2005), one traces light rays backwards from the ob- 
server's image plane to the star's surface, by numeri- 
cally integrating the geodesic equations. The flux emit- 
ted from the point where a light ray meets the star will 
then determine the intensity of that ray. This can be 
called an "observer-centered" approach. For our pur- 
poses, however, it is more convenient to take a "grid- 
centered" approach, since our MHD simulations give us 
the hot spot data on a grid, and it is more convenient 
to calculate the contribution of each grid element to the 
observed flux, and then sum over all grid elements, as 
stated above. So we use the following approach (see, 
e.g., Beloborodov 2002). 

Because of light bending, the light from point B in 
Fig. |2K> needs to be emitted at an angle a such that 
after bending, it travels towards the observer, that is, 
it travels at an angle ip with respect to the original ra- 
dial direction. To calculate the observed flux we need 
a relation between ip and a. We use the Schwarzschild 
metric. From the geodesic equations, we then have (see, 
e.g., Misner, Thorne & Wheeler 1973; Pechenick, Ftaclas 
& Cohen 1983) 

if, - f°° — [1 - — (l _ ^-)} 1/2 (2) 
J R r 2 b 2 r 2 \ r / 

where r is the radial Schwarzschild coordinate, R and 
r g are the star's radius and Schwarzschild radius re- 
spectively, and b is the impact parameter. Using the 
four-velocity of a photon at the point of emission, we 
can relate b to a to get (see, e.g., Beloborodov 2002) 
b = R(l — r g / R)^ 1 / 2 sin a. We now have the desired re- 
lation between ip and a. The integration in equation 
J2J cannot be done analytically, but is relatively easy to 
do numerically. However, the problem in our case is a 
more difficult one, because in the grid-centered approach 
we know ip for each grid element, and we have to solve 
equation JzJ for a. Computationally, this is very difficult 
to do exactly. So we use the cosine relation, an approxi- 
mate relation (due to Beloborodov 2002), which has the 
advantages of being very simple to use and highly accu- 
rate: 

1 -cosa w (1 -cos^O (l- ^) ■ (3) 

In that case the observed flux is given by (Beloborodov 
2002) 

J=(l-^-) 2 /" dS I B (R,a)cosa. (4) 

So the angle ip in equation is replaced by a, and 
we get a prefactor of (1 — r g /R) 2 due to gravitational 
redshift. 

An interesting consequence of light bending is that the 
observer can see some radiation from the far side of the 
star (see, e.g., Beloborodov 2002). In particular, both the 
antipodal hot spots of pulsars can be seen simultaneously 
in some cases. 

3.2. Doppler Effect 

The above discussion is valid if Ie (R, ct) is the inten- 
sity in a reference frame which is at rest with respect 
to the observer. However, our MHD simulations calcu- 
late the intensity in a reference frame which is rotating 



with the star. We need to relate the intensities in these 
two frames. We use the invariance of I^/v^ along a ray 
of light, where I v is the specific intensity, or the inten- 
sity per unit frequency range. Then for a ray emitted 
from a point on the star which is moving with a veloc- 
ity v = (3c and Lorentz factor 7 = (1 — /3 2 ) -1 / 2 , we 
have 1 l v = I' v , /^{l — /3^) 3 , since v and v' are related by 
the Doppler formula v' = uj(l — f3fi). Here /i = cos 9, 
where 9 is the angle, as measured in the unprimed frame, 
between the direction of emission of the light and the 
direction of motion of the emitting surface element of 
the star. The direction of emission of light is given by 
equation ©. Integrating over frequency, we then have 
I = I' /3/i) 4 . Thus in equation ffl for the observed 
flux, we pick up an extra factor of 1/7 4 (1 — Pfi) 4 - Also, 
the areas of the hot spots as measured by photon beams 
in the two frames are related by dS' = 7(1 — f3fi)dS (Ter- 
rell 1959). This is a consequence of the projected area 
dS cos a being a Lorentz invariant. So the flux is now 
given by 

J= { l -T) 2 I dS ' sn 1 R \6 4(R>')cos a. 
v R J Jco Sa >o 7 5 (1-/W 5 

(5) 

3.3. Light Travel Time Effects 

Now we shall take into account the fact that light from 
different parts of the star takes different amounts of time 
to reach the observer. We need to use the general rel- 
ativistic expression for the light travel time, since we 
are dealing with a compact object. We again use the 
Schwarzschild metric. Using the geodesic equations, the 
time difference between light emitted from points A and 
B in Fig. |3Jd reaching the observer is (assuming D 3> R) 




(6) 

We then need to find the apparent position of each grid 
element at a given time. At the time when the point A in 
Fig. |2t> will appear to be at the center of the observer's 
image, the point B will appear to be not where it is shown 
in the figure, but where it was a time At ago, where At 
must satisfy 

St(b At ) = At. (7) 

Here b&t is the impact parameter for light emitted from 
point B a time At ago. This determines the apparent 
position of grid element B at the time when grid element 
A is at the position shown in Fig. EJd. 

In the observer-centered approach, one needs to solve 
equation ©, which is relatively easy to do numerically. 
However in the grid-centered approach, one has to solve 
equation J7J) which, after substituting for (5i(&At) from 
equation ©, becomes an integral equation which can- 
not be solved analytically, and is very difficult to solve 
numerically. So, following Beloborodov(2002), we ex- 
pand the integrand in equation (JSJl in powers of x = 
(1 — cos?/') = (1 — cosa)/(l — r g /R) and then perform 

1 Notation in this subsection: Unprimed quantities arc those 
measured in a frame at rest (with respect to the observer) at the 
star's surface. Primed quantities are those measured in a frame 
rotating with the star. 
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the integration. We keep as many terms as are needed 
for sufficient accuracy (accuracy can be checked by com- 
paring the values obtained from the series with those ob- 
tained by integrating eq. JBJ numerically). The resulting 
equation, although still implicit in At, is much easier to 
solve numerically. We can now calculate the apparent 
position of any grid point at any time, and then com- 
pute the flux as given by equation JSJ by summing over 
all grid elements. 

3.4. Frame Dragging Effects 

The neutron stars considered in our simulations have 
R = f km and M = 1.4M Q . The fastest rotators we 
considered have P = 3 ms. For these parameters, we 
tried to find out the significance of frame dragging by 
numerically evaluating frame dragging corrections to the 
path of light in the equatorial plane of the star, using 
the Kerr metric. We found that corrections to the an- 
gle a (Fig. 03) are at most ~ 4°. Corrections to the 
time delay St (eq. EJ are of the order of a few percent, 
which is not large considering the fact that the effect of 
time delay on the variability curve is itself quite small, 
as we shall see in section 5. The gravitational redshift 
factor has corrections of < 1%. We thus expect the er- 
rors introduced by ignoring frame dragging to be much 
smaller than those introduced by the assumptions in our 
variability model (which we discuss in section 5). This, 
together with the fact that frame dragging effects are 
very difficult to take into account numerically, is why we 
do not do so. Variability curves have been calculated by 
earlier authors taking these effects into account (Braje, 
Romani & Rauch 2000), and they found that the curves 
do not change significantly due to these effects. 

4. EXAMPLES WITH SIMPLE SYNTHETIC SPOTS 
4.1. Point Spot 

It is clear from the foregoing discussion that general 
relativistic effects become stronger with increasing r g , 
and Dopplcr and time delay effects become stronger with 
increasing rotation speed of the star. To see clearly the 
effect of all these effects on the light curve, it is useful 
to consider the simple case of a point spot shown in Fig. 
01 We have only one spot, which is on the rotational 
equator of the star, and the observer is in the equatorial 
plane of the star. When no effects are included, the light 
curve is simply a half-sinusoid (due to the cos ip factor in 
cq. ^) with a maximum when the spot is at point C, and 
minima when the spot is on the far side of the star (Fig. 
|5j) ■ The points labeling the curve correspond to the spot 
being at the respective points in Fig. 01 Let us now look 
at all the effects separately. 

Time delay. Referring to Fig. the light from, say, 
point B, takes longer to reach the observer than that 
from point C. So in the light curve, B shifts to B', and so 
on. Due to this, the light curve gets distorted as shown. 
The peak position of the time-delay curve is arbitrary, 
because the choice of the reference point for the time 
delay (point C in this case) is arbitrary. So time delay 
"bends peaks to the left" . Note that in any case, no shift 
in peak position is observable, since the observer sees 
only one light curve, the one including all effects. We 
discuss the shifts here merely to understand our results 
the better. 
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Fig. 4. — Geometry for the case of a single point spot on the 
equator with the observer in the equatorial plane. 
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Fig. 5. — Schematic light curve for the case shown in Fig. Hlwith 
and without including time delay, over a period of 1 rotation of the 
star. The flux has arbitrary units. 



Doppler. In the above case, the observed flux from 
a certain point is determined only by the costp factor in 
equation Q.When we include the Doppler effect, rela- 
tivistic beaming comes into the picture. In our simple 
case, beaming of radiation towards the observer is max- 
imum at point E, while cos"0 is maximum at point C. 
So the peak occurs at some intermediate point D. The 
Doppler effect thus shifts the peak of the light curve. The 
intensity of the peak also increases dramatically, because 
the spot has a significant velocity along the observer's 
line of sight at that point. For other geometries where 
the spot always moves almost perpendicular to the line 
of sight, beaming reduces the peak intensity. 

General relativity. Without general relativistic ef- 
fects, the visible portion of the star is determined by 
cos-0 > 0. When we include general relativistic effects, 
this condition changes to cos a > 0, with a given by 
equation Since cos a > costp, this condition implies 
that more than half of the stellar surface is visible to the 
observer at any time. For certain geometries and with 
antipodal hot spots, this makes both the spots visible 
simultaneously. In such a case, if the spots are identi- 
cal, the observed flux is almost constant for the dura- 
tion for which both spots are visible (see Beloborodov 
(2002) for a detailed discussion). Even when only one 
spot is visible, the reduced modulation of the flux (due 
to cos a > cosip) flattens the light curve. The peak po- 
sition does not change, however, since the peak occurs 
when cos a is maximum, which is when cos tp is maxi- 
mum. Gravitational redshift decreases the total observed 
flux, except when both the antipodal spots are simulta- 
neously visible, in which case the observed flux can in- 
crease. 
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Fig. 6. — Examples of relativistic and light-travel-time effects on 
the light curves for gaussian spots. Dotted curves do not include 
any effects, dashed curves include only GR effects, dash-dotted 
curves include GR and Doppler effects, and solid curves include all 
effects. The fluxes are normalized. 
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Fig. 7. — Examples of relativistic and light-travel-time effects on 
the light curves of 3-ms pulsars for realistic spots in the PW case. 
The line patterns mean the same as in Fig. [S] 



4.2. Gaussian Spots 

To illustrate the above effects, we show some variability 
curves for spots with a gaussian flux distribution, max- 
imum at the center and tapering outwards. We con- 
sider a hypothetical neutron star with R = 10 km and 
M = 1.4A/©. To bring out clearly the effects of rotation, 
we choose a very small rotation period P = 0.5 ms. The 
spots have a width of 10°. Fig. shows some variability 
curves for different misalignment angles and observer 
inclination angles i. (The inclination angle is the angle 
between the observer's direction and the star's rotation 
axis.) We show the curves over one rotation period of the 
star since, because of the assumptions in our variability 
model, the curves are perfectly periodic with a period 
equal to the star's spin period. Without including any 
effects, the light curve has one or two peaks depending 
on whether one or both the hot spots are seen during 
one rotation period (see R04 for a detailed discussion). 
General relativistic effects reduce the pulse fraction and, 
in the cases where both hot spots are simultaneously vis- 
ible for some time, completely flatten the light curve for 
that duration. The time delay effect is seen to distort the 
light curves. We also see that the Doppler and time de- 
lay effects increase the pulse fraction. These two effects 
are strongest when the hot spots are moving almost di- 
rectly towards or away from the observer, which happens 
at large inclination angles. 

5. VARIABILITY CURVES FOR REALISTIC SPOTS 

We now turn our attention to the realistic spots dis- 
cussed in section 2. To calculate the variability curves, 
we make the following assumptions: only radiation from 
the hot spots contributes to the light curve; when the 
infalling matter falls on the star, its entire kinetic and 
thermal energy is converted into radiation, which is emit- 
ted isotropically; the emitted radiation does not interact 
with the accreting matter. Also, recall that the hot spots 
keep changing, but only slightly, with time. So to model 



the lightcurves, we choose the hot spots at a certain time, 
and then assume those hot spots to be unchanging. Fig. 
\7\ shows the variability curves for the realistic spots in 
the PW case. The neutron star has P = 3 ms. We see 
that general relativistic effects are the strongest. Even 
at such a large rotation speed, the effects of rotation are 
relatively small, which justifies neglecting frame dragging 
effects. 

For clarity, Fig. |H| shows the final light curves (solid 
lines) after including all effects, for different misalign- 
ment and inclination angles. The curves are almost si- 
nusoidal for small inclination angles when Doppler and 
time-delay distortions are relatively weaker, and deviate 
noticeably from a sinusoidal shape for larger inclination 
angles. The pulse fractions depend on the misalignment 
and inclination angles and the shape of the hot spots. In 
most cases the pulse fractions are seen to be quite large 
compared to the observed values of a few percent for real 
stars. We see small pulse fractions in our simulations 
when the misalignment angle or the inclination angle is 
of the order of a few degrees. The most probable reason 
then for the small pulse fractions of real pulsars is that 
they have misalignment angles of the order of a few de- 
grees. Another possible explanation is that scattering of 
the hot spot radiation by the surrounding matter reduces 
the amplitude of oscillations (Brainerd & Lamb 1987). 

We compared the lightcurves from realistic spots with 
those from gaussian ones (dashed lines in Fig. |HJ. We 
noted that the hot spots for the = 30° case are the 
most amenable to approximation by gaussian spots. So 
we approximated them with gaussian spots of the same 
width (34° in this case). The hot spots for the 6 = 60° 
and = 90° cases depart significantly from a round 
shape, but we tried approximating them with gaussian 
spots of the same width as for the = 30° case. The 
gaussian spots are centered at the magnetic poles in each 
case. Their intensity in each case was chosen such that 
the total flux from the star as seen by a distant observer 
looking at the spot from directly above the magnetic pole 
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Fig. 8. — Final light curves, including all effects (solid lines), 
with best-fit gaussian-spot lightcurves (dashed lines), for realistic 
spots in the PW case for 3-ms pulsars. The pulse fractions are 
2% and 8% respectively for the top two panels, ~ 25-30 % for the 
middle two, and ~ 35-40 % for the bottom two. 



would be the same as that in the case of the realistic 
spots. We also set the two antipodal hot spots to be 
identical. The lightcurves thus obtained are also shown 
in Fig. [5] The fractional r.m.s. errors are < 5% for small 
(< 60°) inclination angles and are ~ 5-10 % for large 
(> 60°) inclination angles. The two main reasons for 
gaussian spots not reproducing the lightcurves perfectly 
are that the shapes of the realistic hot spots are more 
complicated than gaussian, and that the two antipodal 
hot spots are not usually identical. Fig. compares the 
variability curves for two identical geometries in the PW 
and non-PW cases. The hot spots, and hence the vari- 
ability curves, have similar shapes in the two cases. We 
found that the flux in the PW case is a few times higher 
than that in the non-PW case in general, as expected 
from the hot spots shown in Fig. [21 However, the pulse 
fractions turned out to be smaller in the PW case. 

We also compared the variability curves of 3-ms and 
5-ms neutron stars for two identical geometries. The 
results are shown in Fig. ^jp. The shapes of the hot 
spots and curves are again similar in these two cases. 
Both the average flux and pulse fraction are higher for 
the 3-ms star. 

We plotted the observed stellar image for our fiducial 
3-ms pulsar in the PW case with and without general 
relativistic effects, for = 90° and i — 15°, shown in 
Fig. ^3 We see that both the antipodal hot spots are 
visible simultaneously in this case. Fig. |3Jd shows that 
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Fig. 9. — (a) Comparison of PW (solid line) and non-PW 
(dashed line) lightcurves, including all effects, (b) Comparison 
of lightcurves for 3-ms (solid line) and 5-ms (dashed line) pulsars, 
including all effects. 



due to bending of light, the star should look larger than 
it actually is, which is also seen here. We do not take the 
Doppler and time delay effects into account in these im- 
ages, since those effects are not noticeable in the images. 

6. CONCLUSIONS 

We modeled the lightcurves of accreting millisecond 
pulsars using hot spots obtained in full 3D MHD simula- 
tions done using the Paczyhski-Wiita potential, for 3 ms 
and 5 ms-period pulsars. We found that the main effect 
of the Paczyhski-Wiita potential is to increase the accre- 
tion rate and the emitted flux. The variability curves in 
our model are strongly affected by general relativistic ef- 
fects, and to a lesser extent by Doppler, time-delay and 
frame dragging effects. General relativistic effects de- 
crease the pulse fraction, while Doppler and light-travel- 
time effects increase it and distort the light curve. The 
amount by which the pulse fraction changes, and the 
distortion of the light curve, depend on the hot spots' 
position and shape, which are determined by the mis- 
alignment angle, and on the observer inclination angle. 
The small pulse fractions observed in real pulsars sug- 
gests that they might have small misalignment angles, of 
the order of a few degrees. 

We compared the lightcurves from the realistic hot 
spots that we obtained from our MHD simulations with 
those from simple hot spots with a gaussian flux distribu- 
tion centered at the magnetic poles. We found that the 
gaussian-spot lightcurves differ from the realistic-spot 
lightcurves by < 10%, and that therefore gaussian spots 
are a reasonable approximation for the realistic ones. We 
plan to investigate in the future how the size and inten- 
sity of these equivalent gaussian spots depend on physical 
parameters pertaining to the star and the disk. 

Numerical variability models, along with analytical 
models, are important tools for studying periodic and 
quasi-periodic oscillations from X-ray pulsars. In the re- 
gion near the star, hot spots, funnel streams and features 
in the accretion disk like density waves could be respon- 
sible for producing these oscillations. Our 3D simula- 
tions are useful for studying these features. However, 
a more accurate model of the variability curves, which 
takes into account emission and absorption of radiation 
by the accreting matter, and temporal changes in the hot 
spots, is needed to determine the role of these features in 
producing the oscillations. Also, one of the major sim- 
plifications in our model is that we ignore the dynamics 
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Fig. 10. — Observed stellar image, without including any effects (left panel) and including general relativistic effects (right panel). Note 
the apparent enlargement of the stellar image and the simultaneous visibility of both antipodal hot spots because of general relativistic 
effects. 



of the matter after it falls on the star's surface. This 
is not a good assumption for millisecond pulsars, where 
thermonuclear burning of the matter falling on the stel- 
lar surface is a possibility (Joss & Li 1980; Bildsten & 
Brown 1997; Bhattacharyya et al. 2005), which could 
change the dynamics of the matter flow and magnetic 
field around the star, especially since the magnetic field 
is relatively weak. In that case, the hot spots would serve 
to give an idea of the place where matter would accumu- 
late and thermonuclear burning could start. 
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